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ABSTRACT 

In  order  to  predict  the  response  time  of  smoke  detectors  in  enclosure  fires,  a 
computational  model  is  developed  for  calculating  the  local  particulate 
concentration  near  the  ceiling.  The  large  scale  smoke  movement  is  approximated 
by  integral  equations  for  the  plume  and  ceiling-jet,  which  originates  in  the 
cold  lower  layer  and  penetrates  into  the  accumulated  smoke  layer  in  the  upper 
portion  of  enclosure.  The  effect  of  coagulation,  which  changes  the  particle 
size  distribution,  is  included  to  enable  predictions  of  an  ionization  smoke 
detector  response.  This  engineering  model  is  designed  to  be  used  in 
combination  with  two- layer  zone  models  for  obtaining  more  detailed  information 
of  smoke  concentration  in  the  upper  layer.  Sample  calculations  have  been  made 
and  comparisons  with  relevant  experimental  data  showed  a reasonable  agreement 
both  in  the  mass  concentration  and  particle  number  concentration  of  smoke  in 
the  ceiling-jet. 


1.  INTRODUCTION 

Comprehensive  fire  models  like  those  of  Mitler  and  Emmons  [1],  and  Jones  [2] 
have  paved  the  way  for  the  deterministic  assessment  of  potential  fire  hazards 
in  buildings.  One  of  the  most  important  elements  in  performing  a hazard 
analysis  is  the  time  of  detection.  Earlier  detection  of  fire  can  provide  a 
better  safety  margin  to  evacuat  building  occupants.  Further,  early  detection 
is  necessary  for  rapid  notification  of  the  fire  department  and  for  rapid 
initiation  of  automatic  fire  suppression.  A great  deal  of  research  has 
already  been  done  to  understand  the  response  of  the  heat  sensitive  devices, 
such  as  automatic  sprinklers,  but  fewer  studies  have  been  made  in  relation  to 
the  response  of  smoke  detectors. 

Two  types  of  smoke  detectors  are  in  common  use  today:  the  photo-electric  type 
and  the  ionization  type.  These  two  types  of  detectors  operate  on  fundamentally 
different  physical  principles.  The  photo-electric  smoke  detectors  use  the 
scattering  of  light  by  smoke  particles.  The  ionization  smoke  detectors  use  the 
reduction  of  ion  current  in  an  ionization  chamber  caused  by  the  attachment  of 
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ions  to  smoke  particles.  Both  theoretical  and  experimental  studies  [3,4]  show 
that  the  response  functions  of  these  detectors  are  strongly  dependent  on  the 
number  concentration  and  the  size  of  smoke  particles.  Thus,  in  order  to 
develop  a predictive  method  for  the  response  of  smoke  detectors,  one  must 
include  the  effects  influencing  the  particulate  component  of  smoke 
characteristics  in  addition  to  a model  for  large  scale  smoke  movement. 

Smoke  aerosol  generated  at  a fire  source  is  carried  upward  as  part  of  a 
buoyant  plume.  The  plume  then  impinges  on  the  ceiling  and  forms  a radial 
outward  flow  under  the  ceiling.  The  nature  of  the  flow  is  generally  turbulent 
with  strong  entrainment;  this  dilutes  the  smoke  aerosol  during  its  way  to  the 
detector.  The  heat  and  particulate  release  rates  are  the  most  important 
factors  influencing  the  smoke  concentration  near  the  ceiling.  Further, 
accumulation  of  the  smoke  layer  in  the  upper  portion  of  the  enclosure  may  vary 
the  effect  of  dilution  thus  affecting  the  detector's  response  time. 

Another  influential  effect  on  the  response  of  smoke  detectors  is  the 
coagulation  of  aerosol.  At  high  concentration,  smoke  particles  frequently 
collide  and  stick  together  as  a result  of  Brownian  motion.  Through  this 
process  the  average  particle  size  increases  while  the  particle  number 
concentration  decreases. 

The  first  attempt  to  model  the  smoke  dynamics  in  an  enclosure  fire  has  been 
made  by  Baum  et  al.  [5].  In  their  model,  the  large  scale  fluid  motion  is 
directly  calculated  from  the  field  equations  of  motion.  Smoke  movement  is 
simulated  by  tracking  mathematical  "blobs"  of  smoke,  and  the  effect  of 
coagulation  is  treated  as  an  internal  process  within  each  blob.  In  a previous 
study,  the  author  has  developed  a model  in  which  local  particle  number 
concentration  is  directly  calculated  from  the  field  equation  of  coagulation 
[6].  However,  there  are  difficulties  in  applying  these  field  equation  models 
to  our  engineering  problems  of  detection,  due  to  their  requirement  for  very 
large  computing  facilities. 

Recently,  Evans  [7]  developed  an  engineering  model  for  the  prediction  of  the 
response  time  of  automatic  sprinklers.  Based  on  the  commonly  used  point  source 
approximation,  known  as  the  Morton-Taylor-Turner  model  [8],  Evans  obtained  an 
analytical  solution  for  the  plume  flow  in  the  two -layer  environment  of 
enclosure  fires.  A noticeable  achievement  in  Evans'  study  is  that  he  showed 
the  prediction  of  response  time  of  ceiling-mounted  detectors  is  possible 
within  the  framework  of  two -layer  zone  models  and  with  sufficient  accuracy. 

Following  Evans'  study,  prediction  of  response  time  of  smoke  detectors  on  the 
basis  of  a zone  model  is  now  possible.  In  this  report,  a method  is  presented 
for  approximating  the  large  scale  smoke  movement  in  plume  and  ceiling-jet, 
which  originates  in  the  cold  lower  layer  and  penetrates  into  the  accumulated 
smoke  layer  in  the  upper  portion  of  the  enclosure.  The  effect  of  coagulation 
is  included  to  enable  predictions  of  the  response  of  ionization  smoke 
detectors.  Results  of  comparison  with  relevant  experimental  data  are  also 
described . 


2 


2.  PLUME  REGION 


2 . 1 Governing  Equations 

Similar  to  Evans'  study,  the  Morton-Taylor-Turner  model  is  used  as  a basis  for 
the  present  plume  model  in  a two- layer  environment. 

The  major  assumptions  of  the  present  model  are: 

(1)  The  fire  plume  is  axisymmetric  and  steady  with  no  swirl; 

(2)  The  environment  is  stagnant  and  stratified; 

(3)  The  density  (temperature)  and  the  particulate  concentration  in  each  layer 
are  uniform; 

(4)  The  vertical  velocity,  temperature  and  particulate  concentration  have 
Gaussian  profiles  in  the  radial  direction; 

(5)  The  ideal  gas  law  is  valid; 

(6)  The  Boussinesq  approximation  is  valid. 

The  profiles  of  velocity,  temperature  and  particulate  concentration  can  be 
written  as  follows: 


u/Ufl,  - exp(-r2/b2),  (1) 

AT/ATm  - ACs/ACsm  - exp  [ - r2/(  Ab) 2 ] . (2) 

where  AT  = T - Ta  and  ACS  = Cs  - Csa  represent  respectively  the  temperature 
difference  and  the  concentration  difference  from  ambient,  and  b(z)  is  the  1/e 
width  of  the  plume  with  respect  to  the  velocity  distribution.  The  profiles  of 
the  temperature  and  the  particulate  concentration  are  assumed  to  be  identical. 

With  the  above  assumptions,  the  governing  equations  for  mass,  momentum,  energy 
and  particulate  concentration  can  be  written  as  follows: 


dz 


r u dr 
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r uz  dr 
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(6) 


where  a is  the  entrainment  constant  determined  experimentally.  Ta  and  Csa 
represent  respectively  the  ambient  temperature  and  the  ambient  concentration 
at  each  height,  and  Tr  is  the  arbitrary  reference  temperature. 

Integration  of  Eqs . (3) -(6)  yields  the  following  relations  governing  plume 

flow  in  a stratified  ambient: 
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(10) 


Eqs.  (7) -(9)  are  equivalent  to  the  equations  used  by  Evans  [7],  and  Eq . (10) 
for  smoke  concentration  has  been  added  in  the  same  form  as  the  energy 
equation.  The  plume  flow  in  the  upper  layer  is  calculated  directly  from  the 
governing  equations  by  imposing  boundary  conditions  at  the  interface  between 
the  upper  and  lower  layer  based  on  plume  flow  in  the  lower  layer.  This  is 
similar  to  one  of  the  several  methods  explored  previously  by  Evans  [7].  In 
deriving  boundary  conditions  at  the  interface,  it  is  assumed  that  the  flux  of 
mass,  momentum,  excess  enthalpy,  and  smoke  are  conserved  across  the  -interface. 
Since  each  layer  may  be  considered  as  having  uniform  temperature  and  smoke 
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concentration  away  from  the  plume  flow  region,  the  right  hand  side  of  Eqs . (9) 

and  (10)  are  non-zero  only  at  the  interface. 

The  plume  flow  in  the  upper  layer  has  to  be  solved  numerically.  In  the  program 
listed  in  APPENDIX,  the  fourth-order  Runge-Kutta  method  is  employed  in 
obtaining  a solution  for  the  plume  flow  in  the  upper  layer. 

2.2  Parameters 

The  values  of  the  entrainment  constant  a and  Gaussian  width  ratio  A are 
necessary  for  evaluation  of  the  plume  model.  From  extensive  experimental 
study,  McCaffrey  [9]  provided  a set  of  correlations  for  the  centerline 
variations  of  temperature  and  velocity  in  buoyant  diffusion  flames. 

McCaffrey's  correlations  for  the  plume  region  (0.20  < z/Q2/-’)  are  consistent 
with  the  analytical  solution  of  Eqs.  (7)  -(9)  for  homogeneous  environment. 

Also,  his  data  are  believed  to  be  the  most  reliable  results  available. 
Therefore,  the  values  of  a and  A were  selected  to  match  the  McCaffrey's 
correlations . 

Comparison  of  McCaffrey's  correlations  for  the  plume  region  with  the 
analytical  solution  of  Eqs.  (7) -(9)  yields  the  following  values  of  a and  A: 


a = 0.118, 

A2  = 1.157.  (11) 


These  values  will  be  used  in  all  calculations  to  appear  in  this  paper.  They 
are  slightly  different  but  very  close  to  the  values  used  by  Zukoski  [10]:  a = 
0.110,  A2  = 1.12.  It  should  be  noted  that  the  vertical  coordinate  z shall  be 
taken  as  measured  from  the  flame  base  to  be  consistent  with  McCaffrey's 
correlations . 

2.3  Comparison  with  Experimental  Data 

Currently,  no  detailed  smoke  concentration  data  for  the  fire  induced  plumes 
are  available.  However,  it  is  meaningful  to  evaluate  this  model  using 
temperature  data,  since  the  smoke  aerosols  are  believed  to  follow  the  same 
mechanism  of  transportation  as  temperature  under  adiabatic  conditions. 

Evans  [7]  provided  one  set  of  experimental  data  for  plume  centerline 
temperatures  in  a two -layer  environment.  Fig.  1 is  a comparison  of  Evans'  data 
with  prediction.  Although  there  is  a sharp  decrease  in  calculated  plume 
temperature  at  the  interface  as  pointed  out  by  Evans,  the  overall  agreement  is 
quite  good.  This  discontinuity  of  temperature  can  be  smoothed  out  by  simply 
assuming  that  the  plume  carries  its  lower-layer  features  to  some  distance  into 
the  upper  layer  (as  shown  by  the  dotted  line  in  Fig.  1). 
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3.  CEILING-JET  REGION 


3.1  Governing  Equations 

Incorporating  the  Boussinesq  treatment  for  buoyancy  and  the  Reynolds'  analogy 
for  convective  heat  transfer,  Alpert  [11]  developed  an  integral  model  for  the 
ceiling-jet  induced  by  a fire  under  an  unconfined  ceiling.  Alpert's 
theoretical  predictions  are  in  good  agreement  with  the  empirical  correlations 
obtained  by  Alpert  himself  [12],  and  Pickard  [13]. 

Assuming  that  the  fire  induced  ceiling- jet  confined  in  an  enclosure  behave  the 
same  way  as  an  unconfined  ceiling- jet  with  the  environment  of  the  upper  layer, 
Alpert's  model  may  be  applied  for  calculating  the  near-ceiling  smoke 
concentration  as  a function  of  radial  distance  from  the  fire  axis. 

The  major  assumptions  are: 

(1)  The  ceiling- jet  is  axisymmetric  and  steady  with  no  swirl; 

(2)  The  environment  is  stagnant  and  uniform; 

(3)  The  velocity,  temperature  and  particulate  concentration  have  half -Gaussian 
profiles  in  the  vertical  direction; 

(4)  The  ideal  gas  law  is  valid; 

(5)  The  Boussinesq  approximation  is  valid. 

A detailed  analysis  of  this  type  of  flow  can  be  found  in  Ref.  [11]..  Thus,  only 
the  final  forms  of  the  governing  equations  and  a brief  description  will  be 
presented  here. 

The  vertical  distribution  functions  of  ceiling- jet  velocity,  temperature  and 
smoke  concentration  can  be  obtained  from  Eq . (1)  and  (2)  by  replacing  the 
coordinate  r by  y,  and  the  length-scale  b by  2 ; H is  the  1/e  thickness  of  the 
ceiling-jet;  implications  for  this  is  that  the  profiles  are  assumed  to  bear 
the  same  relationship  as  in  the  plume.  The  following  are  the  integrated  form 
of  the  governing  equations  for  mass,  momentum,  energy  and  particulate 
concentration,  respectively,  for  an  axisymmetric  ceiling-jet: 
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The  quantities  v,  AT  and  ACS  appearing  in  Eqs . (12) -(15)  are  the  spatially- 
averaged  values  of  the  corresponding  quantities  in  the  ceiling-jet.  The 
relations  between  the  averaged  and  the  maximum  values  are  as  follows: 
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The  quantity  S is  the  shape  parameter  which  depends  on  the  velocity  and  the 
temperature  profiles.  In  the  present  case,  S takes  the  value: 


2 A2 

S = ■ — - 1.0.  (19) 

7T  Ja2/(1  + A2) 
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By  use  of  the  Richardson  number 


R-i 


g h AT 


v" 


(20) 


the  conservation  equations  of  mass,  momentum  and  energy  can  be  combined  into 
two  separate  equations  for  ceiling- jet  thickness  h and  Richardson  number  Rj_ , 
which  results  as  follows: 
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where  a = AT/Ta. 


(22) 


The  entrainment  function  E can  be  expressed  as  a function  of  the  Richardson 
number  R^  as : 


E = Eq  exp[£  (Rio  - R±)] 


(23) 


where  /3  is  an  empirical  constant;  EQ  and  R^Q,  respectively,  are  the  values  of 
E and  R^  at  the  start  of  the  ceiling-jet. 

The  convective  heat  loss  term  Q /(va)  can  be  correlated  to  the  friction  factor 
fw  by  use  of  the  Reynolds'  analogy  between,  the  skin  friction  and  the 
convective  heat  transfer: 
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Kw  = (T  - TW)/(T  - Ta)  characterizes  the  thermal  conditions  of  the  ceiling.  Kw 
= 0 can  be  used  for  the  adiabatic  ceiling,  and  Kw  = 1 can  be  used  for  the 
isothermal  ceiling  at  Tw  = Ta.  The  value  of  Kw  may  be  calculated  using  the 
results  of  other  fire  models  on  wall  and  upper-layer  temperatures. 


Now,  making  use  of  Eqs . (23) -(25),  and  with  proper  boundary  conditions,  Eqs . 
(21)  and  (22)  can  be  solved  numerically.  The  other  quantities  of  interest  can 
be  calculated  by  the  following  relationship : 
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where  the  subscript  o denotes  that  they  are  the  values  at  the  exit  of  the 
stagnation  region,  i.e.  up-stream  boundary  conditions  for  the  ceiling-jet. 

3.2  Parameters 


The  parameters  necessary  for  calculations  in  ceiling- jets  are  the  friction 
factor  fw,  and  EQ  and  /3  used  in  the  entrainment  function.  In  order  to  obtain 
the  local  values  of  fw,  Alpert  calculated  the  thickness  of  the  viscous 
sublayer  formed  between  the  ceiling  and  the  main  portion  of  ceiling-jet. 
However,  as  pointed  out  by  Alpert  himself,  the  value  of  fw  has  little  effect 
on  the  flow  properties.  In  more  recent  experimental  studies  on  fire  induced 
ceiling- jets  [14,15],  values  of  fw  from  0.02  to  0.04  have  been  found  to  match 
the  ceiling  heat- transfer  data.  Further,  You  [16]  pointed  out  that  if  the 
stagnation  point  heat  flux  is  correlated,  then  fw  = 0.025  gives  a best  fit  to 
their  experimental  data  on  ceiling  heat  transfer.  In  this  study,  fw  is  assumed 
constant  and  the  value  0.025  will  be  used. 

Alpert  assumed  that  the  magnitude  of  entrainment  at  the  start  of  the  ceiling- 
jet  is  equal  to  the  plume  entrainment  constant,  and  used  EQ  = 0.12  and  /3  = 

3.9.  If  the  same  assumption  is  applied  to  the  present  model,  then  EQ  becomes 
0.167  (0.118  for  equivalent  Gaussian  profile).  The  value  of  EQ  has  a 
significant  effect  on  the  radial  variation  of  smoke  concentration  and 
temperature,  especially  in  the  distance  range  of  one  or  two  ceiling-heights. 
Further,  the  initial  entrainment  need  not  necessarily  be  equal  to  the  plume 
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entrainment  constant.  Therefore,  several  other  values  are  tried  in  this  study 
and  their  effect  will  be  discussed. 

3.3  Computation 

Alpert  provided  the  upstream  boundary  conditions  for  the  ceiling- jet  in 
relation  to  the  impinging  plume  characteristics,  which  can  be  written  with  the 
present  notations  as. follows: 


ho  = bp/J6> 
vo  = ump/ '■I  ^ > 


ATo  9 ^mp  * 

1 + A2 


^so  0 ^smp 

1 + A2 


sj0  = h (i  + Je/s^)-1 


(29) 


(30) 


The  forth- order  Runge-Kutta  method  is  also  used  for  numerical  computations  in 
the  ceiling-jet. 

3.4  Comparison  with  Experimental  Data 

In  order  to  evaluate  the  ceiling- jet  model,  calculations  for  homogeneous 
environments  are  carried  out  and  the  results  are  compared  with  relevant 
experimental  data  on  unconfined  ceiling- j ets . 

Although  minor  differences  exist  due  to  variations  in  the  plume  related 
constants,  almost  identical  results  to  Alpert 's  calculations  are  obtained  for 
the  ceiling- jet  thickness,  temperature  and  velocity,  on  the  condition  that  the 
same  values  are  used  in  E0  and  f)  (E0  = 0.12,  /3  = 3.9). 

Veldman  et  al . [14]  provided  adiabatic  ceiling- jet  temperature  data  in  their 

experimental  study  on  ceiling  heat  transfer.  Fig.  2 shows  a comparison  of 
their  data  with  prediction.  With  EQ  = 0.12,  ft  = 3.9  and  Kw  = 0 (adiabatic 
condition) , the  calculation  provides  a correlation  which  fits  the  data  quite 
nicely . 
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The  data  of  Veldman  et  al . were  obtained  from  laboratory  scale  experiments; 
heat  release  rates  are  between  1.17  kW  and  1.53  kW  and  ceiling-heights  are 
between  0.584  m and  0.813  m.  On  the  other  hand,  Alpert's  experimental 
correlations  [12]  were  obtained  from  large  scale  tests  with  heat  release  rates 
ranging  from  670  kW  to  100  MW  and  ceiling  heights  from  4.6  m to  13.7  m.  Since 
the  present  study  is  concerned  with  the  very  early  stage  of  an  enclosure  fire, 
it  is  of  great  interest  that  prediction  is  compared  with  full  scale 
experiments  of  small  source  strength  (in  the  range  of  100  kW) . 

Heskestad  [17]  provided  experimental  correlations  obtained  from  fire  tests 
with  heat  release  rates  in  the  range  of  5 kW  to  375  kW  and  the  ceiling  heights 
from  1.2  m to  2.4  m.  Shown  in  Fig.  3 is  a comparison  between  prediction  of 
this  model  with  Heskestad' s correlation  in  radial  variation  of  maximum  excess  > 
temperatures.  The  agreement  is  again  quite  good.  However,  in  obtaining  this 
result,  the  initial  magnitude  of  entrainment  EQ  was  set  to  0.05  which  is  about 
half  of  the  value  used  by  Alpert.  Fig.  4 shows  in  the  same  way  a comparison 
between  the  prediction  and  Heskestad' s correlation  in  maximum  velocities.  The 
experimental  and  computational  conditions  are  the  same  as  in  Fig.  3.  Although 
general  agreement  can  be  seen,  a consistent  over-prediction  is  observed  in 
this  case. 

The  cause  of  inaccuracy  in  the  velocity  variations  is  not  now  understood,  but 
it  may  be  associated  with  the  assumed  Gaussian  velocity  profile.  For  the 
distance  range  greater  than  one  ceiling  height,  the  experimental  correlation 
is  very  close  to  the  calculated  velocity  values  of  equivalent  top-hat 
profiles,  which  is  1 / J 2 times  less  than  the  maximum  velocities,  suggesting 
decay  of  the  Gaussian  profile.  This  inaccuracy  in  velocity  is  not  considered 
important  in  calculating  the  smoke  detector's  response  time,  since  its 
dependence  on  the  gas  velocity  is  generally  quite  small.  Thus,  the  values  EQ  = 
0.05  and  (3  = 3.9  will  be  used  in  the  following  calculations  of  ceiling-jets  in 
room  fire  configurations. 


4.  FILLING  PROCESS  AND  COAGULATION  OF  SMOKE 

At  high  concentration,  the  particle  number  concentration  of  smoke  changes 
quite  rapidly  by  frequent  particle  collisions  and  sticking  together.  This 
effect  of  coagulation  is  specially  important  in  quantifying  the  response  of 
ionization  smoke  detectors,  because  their  response  functions  are  closely 
related  to  the  number  concentration  of  smoke  aerosol.  In  the  room 
configurations,  this  effect  of  coagulation  may  be  estimated  separately  for  the 
three  regions  of  concern:  plume,  ceiling- jet  and  upper  smoke -layer. 
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4.1  Coagulation  in  Plume  and  Ceiling- jet 


Baum  and  Mulholland  [18]  presented  a theory  for  the  coagulation  in  a plume. 
According  to  their  analysis,  the  problem  of  coagulation  in  the  plume  can  be 
transformed  into  a problem  of  homogeneous  coagulation  by  introducing  the 
characteristic  time  Tp  defined  by 


dfp  A2  + 1 

_ = [b (z)  UhjCz)  ] “2  , (31) 

dz  2 


where  b(z)  and  u^z)  , respectively  are  the  Gaussian  width  and  the  maximum 
velocity  in  the  plume  as  defined  in  Eqs . (1)  and  (2)  of  this  paper.  The 
particle  number  flux  0 as  a function  of  height  can  then  be  given  by 

<f>  / 4>0  - 1 / [1  + r a0  r]  , (32) 


where  T is  the  coagulation  coefficient  and  <t>Q  is  the  initial  number  flux. 

By  applying  the  same  analysis  to  the  ceiling-jet,  the  following  expression  can 
be  obtained  for  the  characteristic  time  of  ceiling-jet: 


dr 


\ 


A2  + 1 


[r  i(r)  vm2(r)j "1 


(33) 


Eq . (32)  also  applies  for  estimating  the  effect  of  coagulation  in  the  ceiling - 
jet. 
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4.2  Filling  Process  and  Coagulation 

The  problem  of  quantifying  the  effect  of  coagulation  in  the  upper  layer  can  be 
simplified  by  employing  the  commonly  used  assumption  of  two-layer  zone  models, 
i.e.  the  properties  in  the  upper  layer  are  uniform.  The  equations  describing 
the  conservation  of  mass  and  number  concentration  of  smoke  in  the  upper  layer 
can  generally  be  expressed  as: 


dMs  _ . 

“ ' = msp  ” mse » 
dt 


(34) 


dNf 


dt 


nsp  “ nse  “ Kc * 


V 


Kc  = 


r n„z  dV, 


(35) 


0 

w 


where  Ms  and  Ns , respectively,  are  the  total  mass  and  total  particle  number  of 
smoke  in  the  upper  layer;  mSp  and  nSp  represent  the  source  terms;  mse  and  nse 
represent  the  loss  of  smoke  through  the  doorway;  ns  is  the  local  number 
concentration  of  smoke;  V is  the  volume  of  the  upper  layer. 


With  the  assumption  of  homogeneous  concentration,  the  coagulation  term  can  be 
simplified  as 


K, 


'V 

r ns2  dV  = r Ns2  / V, 


0 


(36) 


and  the  loss  terms  can  be  correlated  to  the  total  mass  M,  and  the  mass  loss 
rate  me  of  the  gas  in  the  upper  layer: 


m 


se 


m. 


M 


M 

1As  > 


n 


se 


m. 


M 


N£ 


(37) 
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Substitution  of  Eqs . (36) -(37)  into  Eqs . (34) -(35)  yields: 


dM, 


dt 


f " 

nu 


= m 


sp 


M 


M 


dN, 


dt 


= n 


sp 


p r 


M 


Nc 


m. 


M 


Nt 


(38) 


(39) 


Expressing  the  equations  in  this  form  has  the  advantages  that  p-P  is 
relatively  insensitive  to  temperature  and  the  model  can  be  easily  combined 
with  other  zone  models  by  which  values  of  M and  me  as  functions  of  time  may  be 
provided.  By  further  substituting  N'  = Ns/M  in  equation  39,  N'  becomes 
insensitive  to  expansion  as  the  gas  is  heated. 

Given  the  values  of  M and  me , Eqs.  (38)  and  (39)  can  be  solved  numerically. 

The  mass  and  number  concentration  of  smoke  in  the  upper  layer  can  then  be 
obtained  by  dividing  the  variables  by  the  upper -layer  volume  V: 


Cs  = Ms  / V, 

ns  = Ns  / V.  (40) 


4.3  Simple  model 


We  conduct  a sample  calculation  here  to  make  a rough  estimate  on  the  effect  of 
coagulation  in  practical  cases  of  early  stage  fire  scenarios.  As  a typical 
case  of  ceiling  height  and  detector  location,  we  choose  H = 3.0  m and  r = 3.0 
m,  and  the  fire  size,  Q = 100  kW.  Eq . (31)  can  be  solved  analytically  for  the 
case  of  uniform  environment,  which  has  the  form: 


f 9 


Tp  = 3 


( + 1)  a 


2/3  f g Q 


>.  Pa.  *3p  > 


-2/3 


(z0-V3  . z ‘ I/3 ) 


(41) 
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Assuming  that  the  flame  height  zQ  is  given  by  the  relation  0.08  = zq/Q2/3 
(McCaffrey  [9]),  and  z by  Eqs . (30),  a value  rp  = 3.3  s2/m3  can  be  obtained 
for  zQ  = 0.5  m,  z = 2 . 84  m and  Q = 100  kW.  Numerical  calculation  of  Eq . (33) 
gives  t- ; = 2.1  s2/m3  at  r = 3 . 0 m for  the  same  source  strength.  The  typical 
value  or  coagulation  coefficient  may  be  T = 1.0  x 10 m3/s  (Lee  and 
Mulholland  [19]).  Now,  substitution  of  the  total  characteristic  time,  r = rp  + 
Tj  into  Eq . (32)  gives  the  ratio  of  initial  and  local  particle  number  flux 

4>  / 4>0,  as  illustrated  in  Table  1. 


Table  1 Effect  of  coagulation  in  plume  and  ceiling-jet:  calculation. 


^o 

(particles/s) 

<t>  / <t>Q 

number  concentration 
at  r - 3 . 0 m 

O 

(particles/mJ) 

r = 0.0  m(r  = rp) 

r = 3 . 0 m 

1.0  x 1012 

0.997 

0.995 

3.57  x 1011 

1.0  x 1013 

0.968 

0.949 

3.41  x 1012 

1.0  x 1014 

0.753 

0.651 

2.34  x 1013 

1.0  x 1015 

0.234 

0.157 

5.64  x 1013 

1.0  x 1016 

0.030 

0.018 

6.46  x 1013 

It  can  be  seen  from  Table  1 that  the  total  effect  of  coagulation  in  plume  and 
ceiling- jet  becomes  quite  important  when  the  particle  number  flux  exceeds  the 
critical  value  of  10 -1-4  particles/sec.  It  also  suggests  that  coagulation  in  the 
plume  is  more  important  than  in  the  ceiling-jet,  because  its  initial 
concentration  is  much  higher  than  in  the  ceiling-jet.  The  effect  of 
coagulation  becomes  even  more  important  when  the  heat  release  rate  has  a 
smaller  value.  The  critical  number  flux  for  a 10  kW  source  is  approximately 
10^3  particles/sec. 


The  coagulation  in  the  upper  layer  also  plays  a significant  role  in  varying 
the  smoke  detectability,  because  it  changes  the  environment  of  plume  and 
ceiling-jet.  If  a time  scale  of  300  sec  is  assumed  for  the  detector's  response 
time,  then  the  coagulation  effect  becomes  important  even  for  the  number 
concentration  of  10^-2  particles/m3;  the  sensitivity  of  the  ionization  detector 
lies  in  this  range. 


It  is  believed  that  the  number  flux  at  the  wall,  where  ceiling- jet  ends,  is 
fed  into  the  upper  layer.  However,  a quantitative  model  for  describing  this 
effect  is  not  currently  available.  Therefore,  more  simplified  model  will  be 
used  here. 
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In  the  simple  model,  the  effect  of  coagulation  in  plume  and  ceiling- jet  is  not 
directly  calculated,  but  included  in  the  source  term  of  number  concentration 
in  the  upper  layer.  If  the  coagulation  in  the  ceiling- jet  is  neglected,  then 
the  effective  particle  number  flux  may  be  calculated  by  Eqs . (31)  and  (32), 
assuming  the  characteristic  time  r does  not  change  significantly  by  the 
presence  of  the  accumulated  upper  smoke  layer.  By  numerically  integrating  Eq . 
(39)  using  the  effective  number  flux,  the  average  number  concentration  in  the 
upper  layer  can  be  obtained  as  a function  of  time.  Now,  using  the  effective 
number  flux  as  a source  term  for  the  plume,  the  local  number  concentration  in 
the  ceiling- jet  can  be  calculated  by  the  same  method  described  in  Sec.  2 and  3 
of  this  paper.  Note  that  if  the  effect  of  coagulation  is  neglected,  the 
equations  for  the  number  concentration  become  identical  with  the  equations  for 
the  mass  concentration  of  smoke.  The  exact  form  of  equations  for  the  number 
concentration  can  be  obtained  by  replacing  the  variable  Cs  by  ns  in  Eqs.  (6), 
(10),  (15)  and  (28). 

Based  on  the  models  described  in  Sec.  2,  3 and  in  this  section,  a subprogram 
has  been  developed  for  calculating  the  temperature,  and  mass  and  number 
concentration  of  smoke  in  ceiling- jet  in  a two- layer  environment.  The  program 
needs  input  of  the  upper- layer  thickness,  average  temperature  and  mass  loss 
rate  of  gas  through  the  doorway,  which  should  be  calculated  by  other  two -layer 
zone  models.  The  current  version  of  the  program  is  designed  to  be  used  in 
combination  with  the  program  FAST.  FAST  is  a comprehensive  multi -compartment 
fire  model  whose  detailed  description  can  be  found  in  Ref.  [2,20].  A complete 
list  of  the  subprogram  can  be  found  in  APPENDIX  of  this  paper. 


5.  COMPARISON  WITH  EXPERIMENTS  IN  ROOM  CONFIGURATIONS 

Hochiki  Corp.  of  Tokyo,  Japan  provided  temperature  and  smoke  concentration 
data  collected  in  their  test  room  for  the  purpose  of  evaluating  the 
performance  of  smoke  detectors.  The  size  of  the  room  is  6.0  m D.  x 10.0  m W.  x 
4.0  m H.  The  effective  ceiling-height  can  be  changed  to  3 . 0 m by  placing 
additional  floor-boards  1.0  m above  the  original  floor.  The  walls  are  made  of 
hollow  concrete  blocks  of  25  cm  thick,  whose  effective  thickness  may  be  about 
5 cm.  The  ceiling  is  made  of  1.6  cm  thick  gypsum  plaster  boards  suspended  from 
the  concrete  floor  slab.  The  doorways  and  exhaust  vents  are  sealed  during  the 
tests  so  that  the  only  effective  opening  is  a 30  cm  x 30  cm  drain  hole  located 
at  the  center  of  the  floor.  17  thermocouples  and  12  extinction  meters  are 
permanently  installed  in  the  test  room  to  measure  the  radial  and  vertical 
distribution  of  temperature  and  smoke  concentration.  Additionally,  one 
Measuring  Ionization  Chamber  (MIC)  [21]  is  mounted  on  the  ceiling  3.0  m from 
the  center.  The  MIC  is  a standard  measuring  device  of  smoke  concentration  used 
for  testing  the  sensitivity  of  ionization  smoke  detectors.  The  MIC  employs  the 
same  physical  principle  as  ionization  smoke  detectors.  The  extinction  meters 
used  for  the  experiments  have  equivalent  sensitivity  to  UL  Standard 
specifications  [22].  The  extinction  coefficient  K can  be  obtained  by  the 
following  relation  from  the  meter  readings: 
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In 


(42) 


K = 


L 


where  L is  the  path  length,  IQ  the  intensity  of  incident  light,  and  I the 
intensity  of  transmitted  light. 

The  most  generally  used  combustion  materials  for  testing  the  performance  of 
smoke  detectors  are:  cotton,  heptane,  polyurethane  and  wood  cribs.  Among 
these,  the  heptane  pool  fire  has  been  selected  for  making  comparisons,  since 
it  can  be  regarded  as  a constant  heat/smoke  source  with  the  exception  of  about 
1 min  of  transient  period  after  ignition. 

The  size  of  the  pan  used  for  the  experiments  is  33  cm  x 33  cm  and  5 cm  deep, 
which  was  positioned  on  the  floor  at  the  center  of  the  room.  From  mass  loss 
data,  the  average  heat  release  rate  has  been  found  to  be  approximately  85  kW. 
The  extinction  coefficient  and  MIC  output  taken  at  the  point  3.0  m from  the 
fire  axis  and  5 cm  below  the  ceiling  is  used  for  comparison  with  prediction  of 
ceiling-jet  calculations. 

Using  the  subprogram  listed  in  APPENDIX  in  combination  with  the  program  FAST, 
calculations  have  been  made  for  the  two  different  ceiling  heights,  3.0  and  4.0 
m,  for  which  experimental  data  were  available.  The  heat  release  rate  was 
assumed  to  rise  exponentially  in  the  first  1 min  after  ignition  and  then 
become  constant  (85  kW)  after  this  transient  period.  The  radiation  loss 
fraction  from  the  flame  zone  was  assumed  to  be  0.30,  which  was  subtracted  from 
the  heat  release  rate  to  obtain  the  effective  heat  flux  of  the  plume.  The  mass 
generation  rate  of  particulate  ihp  was  assumed  proportional  to  the  mass  flux  of 
fuel,  and  the  value  of  0.0130  (Mulholland  [23])  was  used  for  the  particle 
conversion  ratio  e = nip/iiif,  where  ihf  is  the  mass  flux  of  fuel.  The  particle 
number  flux  was  calculated  by  assuming  a volume  equivalent  diameter  Dv  of  0.16 
^m  and  the  particulate  density  ps  of  2.0  x 103  kg/m 3: 


n 


sp 


= msp  / (1/6  * 7r  Dv3  ps)  . 


(43) 


The  average  particle  size  Dv  has  been  determined  by  the  experimental  values  of 
extinction  coefficient  and  MIC  output.  The  exact  relation  for  obtaining  the 
volume  equivalent  diameter  will  be  shown  in  the  latter  part  of  this  section. 
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There  is  a well  accepted  relation  between  the  mass  concentration  of  smoke  and 
the  extinction  coefficient  during  flaming  combustion: 


K = 8 Cs , 


(44) 


where  5 is  a dimensional  coefficient  equal  to  7.6  x 10^  m^/kg,  which  was  found 
out  experimentally  by  Seader  and  Einhorn  [24],  and  Lee  and  Mulholland  [19]. 
This  relation  is  used  to  interpret  the  calculated  mass  concentration  into  the 
extinction  coefficient. 

The  output  of  MIC  is  usually  expressed  by  the  quantity  called  Y-value,  which 
was  introduced  by  Hosemann  [3]  in  his  theoretical  study  on  ionization 
chambers : 


Y = 


i 


(45) 


where  iQ  is  the  initial  ionization  current  across  the  positive  and  negative 
electrodes  and  i is  the  ionization  current  in  presence  of  smoke  particles. 
According  to  his  analysis,  the  Y-value  is  proportional  to  the  product  of  the 
particle  number  concentration  and  effective  particle  diameter  ns • d as: 

ns  d = rj  Y,  (46) 


where  rj  is  the  chamber  constant  defined  by 


J<*i  q 

V = 3 — . (47) 

CB 

In  Eq.  (47),  q is  the  ion  generation  rate  by  the  radio-active  source,  a^_  is 
the  recombination  coefficient  of  ions  and  Cg  is  the  Bricard  attachment 
coefficient  equal  to  0.3  cm^/sec . The  value  of  I/77  characterizes  the 
sensitivity  of  an  ionization  chamber.  In  an  experimental  analysis,  Helsper  et 
al . [25]  showed  that  the  value  of  1/rj  for  a Measuring  Ionization  Chamber  is 

approximately  3.5  x 10" ^ m^ . This  value  is  used  for  making  comparisons  between 
the  calculated  number  concentration  and  the  MIC  output. 
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Eqs , (44) -(47)  have  also  been  used  to  calculate  the  average  particle  size  Dv 
of  heptane  smoke  from  the  experimental  values  of  extinction  coefficient  and 
MIC  output.  Given  the  instantaneous  values  of  extinction  coefficient  and  Y- 
value , the  volume  equivalent  diameter  Dv  can  be  obtained  by  the  following 
relation: 


Dv  " 


K 


Y <5  *7  ( 1/6  • 7r  ps)  " 


1/2 


(48) 


In  obtaining  the  relation  of  Eq . (48),  it  has  been  assumed  that  .the  volume 

equivalent  diameter  Dv  is  equal  to  the  effective  average  diameter  d used  in 
the  response  function  of  ionization  smoke  detectors.  The  only  unknown, 
parameter  in  Eq . (48)  is  the  ratio  of  extinction  coefficient  to  Y-value,  K/Y. 
This  parameter  can  be  understood  as  a characteristic  quantity  influenced  by 
the  size  and  shape  of  smoke  particles.  The  exact  value  of  K/Y  depends  on  the 
materials  and  the  conditions  of  combustion. 

In  order  to  verify  the  calculations  by  FAST,  the  calculated  upper- layer 
temperatures  are  compared  with  the  temperature  data  taken  at  the  points  near 
the  corners  of  the  room  and  5 cm  below  the  ceiling.  Fig.  5 and  6 show 
comparisons  of  calculated  and  measured  upper- layer  temperatures  for  two 
different  ceiling  heights.  The  experimental  data  shown  are  the  average  of 
temperatures  taken  at  four  different  points;  the  distance  of  each  point  from 
the  fire  axis  is  identical  (4.5  m)  and  its  distance  from  the  nearest  wall  is 
approximately  1.0  m.  For  both  cases  of  ceiling  height,  the  agreement  between 
the  prediction  and  measurement  is  quite  good.  From  the  fact  that  the 
temperatures  have  been  measured  near  the  ceiling,  this  may  explain  why  the 
upper-layer  temperatures  are  somewhat  over-predicted.  However,  it  may  not  pose 
an  important  effect  on  the  calculations  of  smoke  concentration. 

Fig.  7 shows,  in  the  same  way,  a comparison  between  the  calculated  and 
measured  mass  concentration  of  smoke  in  terms  of  extinction  coefficient.  The 
calculated  result  shown  here  is  the  output  of  FAST,  thus  representing  the 
average  mass  concentration  of  smoke  in  the  upper  layer.  The  experimental  data 
are  for  the  point  5 cm  below  the  ceiling  and  3.0  m from  the  fire  axis. 

Although  an  qualitative  agreement  can  be  seen  between  the  two,  there  is  a 
significant  difference  in  a quantitative  sense.  This  suggests  that  the  smoke 
concentration  in  the  ceiling- jet  is  much  higher  than  the  averaged  values  in 
the  upper  layer.  Similar  results  have  been  obtained  also  in  the  case  of  4.0  m 
ceiling  height  (Fig.  8). 

Fig.  9 shows  comparisons  of  ceiling-jet  calculations  with  the  experimental 
data  for  the  mass  concentration  of  smoke.  The  experimental  data  shown  here  are 
of  the  same  location  as  in  Fig.  5,  and  the  calculated  results  are  for  the  same 
radial  position  in  the  ceiling-jet.  Although  there  are  still  small  differences 
between  the  predictions  and  experimental  data,  the  agreement  is  quite 
reasonable  taking  into  account  the  fact  that  there  is  a 10  to  20  % uncertainty 
in  the  validity  of  Eq . (44)  [19]. 
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Finally,  Fig.  10  shows  a comparison  of  calculated  and  measured  number 
concentration  of  smoke  in  the  ceiling- jet  in  terms  of  MIC  output.  Similar  to 
the  results  in  Fig.  6,  the  calculation  provides  somewhat  lower  values  than  the 
experimental  data.  The  general  agreement  is  again  quite  reasonable.  The 
decreasing  rate  of  change  of  Y with  respect  to  time  is  a result  of  a 
decreasing  number  concentration  resulting  from  particle  coagulation.  It  should 
be  noted  however  that  the  assumption  to  set  Dv  and  d equal  may  not  be 
appropriate,  since  the  exact  relation  between  the  two  size  parameters  is  not 
known  for  the  smoke  aerosol  from  flaming  combustion,  which  usually  has  quite  a 
complicated  agglomerated  shape.  This  leaves  an  uncertainty  in  the  absolute 
value  of  calculated  number  concentration. 


CONCLUSION 

A method  is  described  for  calculating  the  local  particulate  concentration  near 
the  ceiling  in  an  enclosure  fire.  The  large  scale  smoke  movement  is 
approximated  by  integral  equations  for  the  plume  and  ceiling-jet,  which 
originates  in  the  cold  layer  and  penetrates  into  the  accumulated  upper  smoke 
layer.  The  effect  of  coagulation  is  treated  by  a simple  two- layer  zone  model. 
The  key  source  parameters  are  the  fuel- to-particulate  conversion  ratio,  the 
volume  equivalent  diameter  of  particulate  and  the  heat  release  rate.  Sample 
calculations  have  been  made  and  comparisons  with  relevant  experimental  data 
showed  a fairly  encouraging  result. 

Direct  measurement  of  number  concentration  in  high  density  smoke,  and  a more 
detailed  analysis  on  the  response  of  smoke  detectors  especially  for  the 
particle  of  agglomerated  shape,  are  desirable  for  validation  and  refinement  of 
the  model. 

Presently,  application  of  the  model  is  limited  to  a fire  source  of  flaming 
combustion.  Refinements  must  be  made  before  applying  this  model  to  smoldering 
combustion,  which  often  precedes  the  flaming  combustion  in  an  early  stage  fire 
scenario . 
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NOMENCLATURE 


a AT/Ta 

b 1/e  width  of  plume  [Eq.  (1)] 

Cg  Bricard  attachment  coefficient 

Cp  specific  heat 

Cs  particulate -mass  concentration 

ACS  particulate-concentration  difference  from  ambient 
d effective  average  diameter  of  particulate  [Eq.  (46)] 

Dv  volume  equivalent  diameter  of  particulate  [Eq.  (48)] 

E entrainment  function  of  ceiling- jet 

fw  ceiling  friction  factor 

g gravitational  acceleration 

h ceiling- jet  thickness  for  equivalent  top-hat  profile 

hw  ceiling  heat- transfer  coefficient 

I light  intensity 

i ionization  current 

K extinction  coefficient  [Eq.  (42)] 

Kc  loss  of  particles  by  coagulation  [Eq.  (36)] 

Kw  parameter  of  thermal  boundary  condition  at  ceiling  [Eq.  (25)] 
L path  length  of  extinction  meter 

Z 1/e  thickness  of  ceiling-jet 

M total  mass  of  gas  in  upper  layer 

m mass  flux  of  gas 

Ms  total  particulate-mass  in  upper  layer 
ms  particulate -mass  flux 

Ns  total  number  of  particles  in  upper  layer 
ns  particulate -number  concentration 

ns  particulate -number  flux 

Pr  Prandtl  number 

Q heat  release  rate 

Q parameter  of  ceiling  heat  transfer  [Eq.  (17)] 

q ion  generation  rate 

R-l  Richardson  number 

r radial  distance  from  fire  axis 

S shape  parameter  of  ceiling-jet 

St  Stanton  number  [Eq.  (24)] 

T absolute  temperature 

AT  temperature  difference  from  ambient 

u vertical  velocity  of  plume 

V volume  of  upper  layer 

v radial  velocity  of  ceiling-jet 

Y Y-value  of  MIC  output  [Eq.  (45)] 

y vertical  distance  measured  from  ceiling 

z vertical  distance  measured  from  plume  source 
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Greek  letters 

a plume  entrainment  constant 

ion  recombination  coefficient 
/3  parameter  of  entrainment  [Eq.  (23)] 

F coagulation  coefficient 

S constant  defined  by  Eq . (44) 

€ fuel- to-particulate  conversion  ratio 

r?  chamber  constant  of  ionization  chamber 

A Gaussian  width  ratio  [Eq.  (2)] 

p density 

r characteristic  time  of  coagulation  [Eq.  (31),  (33)] 

rw  ceiling  shear  stress 

<j>  particle  number  flux 

Superscripts 

_ spatially  averaged  quantity 

Subscripts 
a ambient 

e mass  flux  through  doorway 

f fuel 

i interface 

j ceiling- j et 

1 lower  layer 

m maximum  value 

o initial  value 

p plume 

r reference  state 

s smoke  particulate 

u upper  layer 

w wall 
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Fig.  1 Comparison  of  experimental  data  with  prediction  in  plume  centerline 
temperatures.  Solid  and  dotted  lines,  calculation;  circles,  Evans  [7]. 


1.17-1.53  kW 
0.584-0.813  m 


H II 


CD 


CM 

O 


3 

E 

•r*l  • 

x i— i 
CO  cO 

e 

w 4J 
0) 


o 

•rH 

e 

x> 

cO 

cQ 

E 

GO 

X! 

cO 

"O 

f— 1 

e 

•rH 

0) 

o 

'O 

CO 

> 

c 

w 

•H 

0) 

rH 

c 

o 

O 

x 

eH 

*rH 

X 

o 

o 

•H 

. o. 

T3 

e 

0)  o 

Xf  'H 

ex  4J 
cO 

-C  i— i 

u 3 

•H  On 

5 r— I -V 
CO  <rt 
n3  O <- 
X>  /■ 


cO 

"d 


“<N 

0)  C 
C 


rX 

cO 

X) 

rH  E- 

V. 

d 

C 

CD 

'O 

•X  *M 

E 

rH 

®r->l 

O m 

X 

00  I 

CD 

ex 

• b 

X 

W CM 

CD 

(D 

X <5 

4-1 

3 03 

O 

X 

cO  U 
C ^ — 
O <D 
w ex  i 
•H  E H 
x o < 

CO  U 

ex  ii 
E x) 
o <d  * 

O *r— ) I 

1 H 

!a0  < 

CM  C 


£>0  *H  <f 
•H  CD  rX 

Cl,  O — ' 


CO 


CD 


CM 


♦^IV 


26 


5-375  kW 
1.2- 2.4  m 


OON.CDin^COCNJ'^-O 


♦ ^IV 


27 


Fig.  3 Comparison  between  prediction  and  experimental  correlation  in  maxim' 
excess  temperatures  of  ceiling- jet.  Solid  line,  calculation;  dashed  line, 
Heskestad  [17].  ATm“  = ATm  [C  2pa2 g H5  / (Ta  Q2 ) ] 1 ' • 
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Fig,  5 Comparison  of  calculated  and  measured  upper- layer  temperatures 
enclosure.  Solid  line,  calculation^  daslied  line,  experiment.  H 3.0m 
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Fig.  7 Comparison  of  calculated  upper- layer  smoke  concentration  with  measured 
ceiling- jet  smoke  concentration.  Solid  line,  calculation;  dashed  line, 
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Fig.  8 Comparison  of  calculated  upper- layer  smoke  concentration  with  measured 
ceiling- jet  smoke  concentration.  Solid  line,  calculation;  dashed  line, 
experiment.  H = 4.0  m. 
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Fig.  9 Comparison  between  calculated  and  measured  ceiling- jet  smoke 
concentration.  Solid  lines,  calculations;  dashed  lines,  experiments. 
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Fig.  10  Comparison  between  calculated  and  measured  particulate  number 
concentration  in  ceiling-jet.  Solid  line,  calculation;  dashed  line, 


APPENDIX  A.  PROGRAM  LISTING 


C 

SUBROUTINE  SMKINIT 
C 

COMMON /SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 

. RDCT , RIO , SP , ST , TO 
C 

REAL  KK , LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL 
C 

ALFA  = 0.118 
GCP  = 1.007 
GRV  = 9.80665 
GRO  = 1.2046 
LMD2  = 1.157 
TO  = 293.0 
C 

BETA  =3.9 
EJO  =0.05 
FW  =0.025 
KK  =1.0 
PR  = 0.7030 
RDCT  =3.0 
SP  =1.0 

ST  = FW  * PR**(-2./3.)  / 2. 

C 

COAG  = 1.0E-15 
PDIA  = 0.16E-06 
RCOAG  = GRO  * COAG 
RSMK  = 2.0E+03 
C 

DLTT  =1.0 
PI  = 3.141592654 

PMAS  = 1./6.  * PI  * (PDIA**3)  * RSMK 
C 

MSTTL  =0.0 
NSTTL  =0.0 
C 

WRITE (7,700) 

C 

RETURN 

700  FORMAT ("  TIME  ",5X,"TU  ",5X,"DTP  ",5XS"DTJ  " ,5X,"CSU  " ,5X, 

. "CSP  " , 5X , "CSJ  " , 5X, "NSU  ",5X,"NSP  ",5X,"NSJ  " ,5X, 

. "UP  " , 5X, "VJ  ",/) 

END 
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c 

SUBROUTINE  SMKTPORT 
C 

C COMMON  BLOCK  AND  INITIALIZATION  FOR  THE  "FAST"  MODEL 

C 

C 

C NR  = MAXIMUM  NUMBER  OF  COMPARTMENTS 

C NN  = NUMBER  OF  NODES  IN  THE  SEPARATORS  (WALLS  AND  CEILINGS) 

C NT  = THE  MAXIMUM  NUMBER  OF  EQUATIONS  TO  BE  SOLVED  (4  * NR) 

C NV  = MAXIMUM  TIME  INTERVALS 

C NS  = MAXIMUM  NUMBER  OF  SPECIES  TO  BE  TRACKED 

C NWALL  = NUMBER  OF  DESCRETE  WALL  SURFACES  (CEILING,  UPPER  WALL  . ..) 

C MXSLB  = MAXIMUM  NUMBER  OF  DIFFERENT  MATERIALS  IN  A WALL 

C 

PARAMETER  (NR=11,  NN  = 24 , NV  = 21,  NS  = 10) 

PARAMETER  (NT  = 4*NR+2*NR*NS , MXSLB=3 , NWAL=4) 

COMMON/F2C17 A/GAMMA , G , SIGM , CP , TA , RA , PREF , RGAS , POFSET , PA , TREF , 

. S S ( NR , NR , 4 ) , SA(NR,NR,4) , AS (NR, NR, 4) , AA ( NR , NR , 4 ) , 

. S AU ( NR , NR , 4 ) ,ASL(NR,NR,4) , NEUTRAL (6 ,NR,NR) , MINMAS , NLSPCT , 

. QF(2 , NR) ,QR(2,NR) ,QC(2,NR) , QFR(2 ,NR) ,QFC(2,NR) , 

. EMP(NR) , EMS (NR) , EME(NR) , APS (NR) ,NWV(NR,NR) , RADIUS (NR , NR) , IVERS , 

. DDT , LFMAX , QEXP , LFBO , LFPOS , LFBT , QRADRL , H COMB A , SWITCH ( 10 ) , 

. BFIRED(NV) , AFIRED(NV) , HFIRED(NV) , TFIRED(NV) , TFMAXT , TFIRET , IFIRED , 

. MPRODR(NV , NS ) , MFIRET(NS) , P(NT) , PMIN(NT) , ZMIN(NR) , 

. BW(NR,NR,4) , HH ( NR , NR , 4 ) , HL ( NR , NR , 4 ) ,NW(NR,NR) ,WINDC(NR) , 

. NOPNMX , NRFLOW , HHP ( NR , NR , 4 ) , HLP(NR, NR, 4) , WINDV , WINDRF , WINDPW , 

. DELTAT , LPRINT , NSMAX , LDIAGP , LDIAGO , ITMMAX , MAXI NR , 

. 1ST , ITMSTP , ID I AG , DERIV(NT) , STIME , NM1 , N , N2 , N3 , N4 , 

. BR(NR) , DR(NR) , HR(NR) , AR(NR) , HRP(NR) ,VR(NR) , HRL(NR) , 

. PW , PM , WC , WH , WO , QP , TE , PPMDV ( 2 , NR , NS ) , TAMB(NR) , RAMB (NR) , PAMB(NR) , 

. FKW (MXSLB , NWAL , NR) , CW (MXSLB , NWAL, NR) , RW(MXSLB , NWAL , NR) , 

. FLW( MXSLB, NWAL, NR) , EPW(NWAL, NR) , NDIV( MXSLB , NWAL, NR) , NSLB(NWAL, NR). 
. , TWJ ( NWAL , NR , NN ) , TWE(NWAL,NR) , QSRADW(2 , NR) ,QSCNV(2,NR) , HFLR(NR) , 

. MASS (2 ,NR,NS) , TOXICT(2 ,NR,NS) , PPM(2 ,NR,NS) , ACTIVS(NS), 

. NCONFG , NDUMPR , LCOPY , CFILE , DFILE , TITLE , TERMXX 
CHARACTER  CFILE(5)*17,  TITLE(50)*1,  DFILE*17 
REAL  MASS , MPRODR , MFIRET , MINMAS , LC50CF , NEUTRAL 
LOGICAL  ACTIVS, SWITCH 
C 

COMMON/SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MGAS , MG OUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 

. RDCT , RIO , SP , ST , TO 
C 

REAL  KK , LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL 
C 

DIMENSION  WRK(2 ) 

EXTERNAL  FSMK 

C 

OLDMAS  = RAMB(l)  * TAMB(l)  / P(N+1)  * P(N2+1) 

C 

MGAS  = MAX ( OLDMAS , MINMAS) 

MGOUT  = SS( 1,2,1)  + SA( 1,2,1) 

MSDT  = MFIRET ( 9 ) 

NSDT  = MFIRET ( 9 ) / PMAS 
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c 

WRK(l)  = MSTTL 
WRK(2)  = NSTTL 
C 

CALL  RKCAL  (FSMK,  TIME,  DLTT , 2,  WRK) 
C 

MSTTL  = WRK(l) 

NSTTL  = WRK(2) 

C 

RETURN 

END 

C 
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SUBROUTINE  FSMK  (F,  X,  Y) 

C 

DIMENSION  F( 2 ) ,Y(2) 

C 

COMMON/SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 
. RDCT , RIO , SP , ST , TO 
C 

REAL  KK , LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL 
C 

F(l)  = MSDT  - Y ( 1 ) * (MGOUT /MGAS ) 

F( 2 ) = NSDT  - RCOAG/MC-AS*(Y(2)**2)  - Y(2)*(MGOUT/MGAS) 

C 

RETURN 

END 

C 
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SUBROUTINE  SMKRESLT  (TIME) 


C 

C COMMON  BLOCK  AND  INITIALIZATION  FOR  THE  "FAST"  MODEL 

C 

C 

C NR  = MAXIMUM  NUMBER  OF  COMPARTMENTS 

C NN  = NUMBER  OF  NODES  IN  THE  SEPARATORS  (WALLS  AND  CEILINGS) 

C NT. = THE  MAXIMUM  NUMBER  OF  EQUATIONS  TO  BE  SOLVED  (4  * NR) 

C NV  = MAXIMUM  TIME  INTERVALS 

C NS  = MAXIMUM  NUMBER  OF  SPECIES  TO  BE  TRACKED 

C NWALL  = NUMBER  OF  DESCRETE  WALL  SURFACES  (CEILING,  UPPER  WALL  ...) 

C MXSLB  = MAXIMUM  NUMBER  OF  DIFFERENT  MATERIALS  IN  A WALL 

C 

PARAMETER  (NR=11,  NN  = 24,  NV  = 21,  NS  = 10) 

PARAMETER  (NT  = 4*NR+2*NR*NS , MXSLB=3 , NWAL=4) 

COMMON/ F2C17A/GAMMA , G , SIGM , CP , TA , RA , PREF , RGAS , POFSET , PA , TREF , 

. S S ( NR , NR , 4 ) , SA(NR , NR , 4) ,AS(NR,NR,4) , AA ( NR , NR , 4 ) , 

. SAU(NR,NR, 4) , ASL(NR , NR , 4) , NEUTRAL (6 , NR , NR) , MINMAS , NLSPCT , 

. QF(2 ,NR) , QR( 2 , NR) ,QC(2,NR) , QFR(2 ,NR) ,QFC(2,NR) , 

. EMP (NR) , EMS (NR) , EME(NR) , APS (NR) ,NWV(NR,NR) , RADIUS (NR, NR) , IVERS , 

. DDT , LFMAX , QEXP , LFBO , LFPOS , LFBT , QRADRL , HCOMBA , SWITCH ( 10 ) , 

. BFIRED(NV) , AFIRED(NV) , HFIRED(NV) , TFIRED(NV) , TFMAXT , TFIRET , I FIRED , 
. MPRODR(NV,NS) , MFIRET(NS) , P(NT) , PMIN(NT) , ZMIN(NR) , 

. BW ( NR , NR , 4 ) , HH ( NR , NR , 4 ) , HL ( NR , NR , 4 ) , NW(NR , NR) ,WINDC(NR) , 

. NOPNMX , NRFLOW , HHP ( NR , NR , 4 ) ,HLP(NR,NR,4) ,WINDV,WINDRF,WINDPW, 

. DELTAT , LPRINT , NSMAX , LDIAGP , LDIAGO , ITMMAX , MAXINR , 

. 1ST , ITMSTP , IDIAG , DERIV(NT) , STIME , NM1 , N , N2 , N3 , N4 , 

. BR(NR) , DR(NR) , HR (NR) , AR(NR) , HRP (NR) ,VR(NR) , HRL(NR) , 

. PW , PM , WC , WH , WO , QP , TE , PPMDV ( 2 , NR , NS ) , TAMB(NR) , RAMB(NR) , PAMB(NR) , 

. FKW (MXSLB , NWAL , NR) , CW(MXSLB , NWAL, NR) , RW (MXSLB , NWAL , NR) , 

. FLW( MXSLB, NWAL, NR) , EPW(NWAL, NR) , NDIV (MXSLB , NWAL, NR) , NSLB(NWAL, NR) 
. , TWJ (NWAL , NR , NN) , TWE(NWAL,NR) , QSRADW( 2 , NR) ,QSCNV(2,NR) , HFLR(NR) , 

. MASS (2 ,NR,NS) , TOXICT(2 ,NR,NS) , PPM(2 ,NR,NS) , ACTIVS(NS), 

. NCONFG , NDUMPR , LCOPY , CFILE , DFILE , TITLE , TERMXX 
CHARACTER  CFILE(5)*17,  TITLE ( 50 )*1,  DFILE*17 
REAL  MASS , MPRODR , MFIRET , MINMAS , LC50CF , NEUTRAL 
LOGICAL  ACTIVS , SWITCH 
C 

COMMON/SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MG AS , MGOUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 

. RDCT , RIO , SP , ST , TO 
COMMON/HFOT/HFOT 
C 

REAL  KK , LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL 
REAL  NSJ , NSP , NSU 
C 


BP 

= 0, 

,0 

CSP 

= 0. 

,0 

DTP 

= 0. 

,0 

NSP 

= 0. 

.0 

UP 

= 0, 

.0 

CSJ 

= 0, 

0 

DTJ 

= 0. 

,0 
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HJ  =0.0 
NSJ  =0.0 
VJ  =0.0 
C 

CSU  = MSTTL  / P(N2+1) 

NSU  = NSTTL  / P(N2+1) 

C 

TL  = P(N3+1) 

TU  = P(N+1) 

TDF  = TU  - TL 

TW  = TWJ (1,1,1) 

QPLM  = QFC (1,1)  + ( -QP+CP*(TE-TL) )*EMP(1) 

ZJO  = ( HR ( 1 ) - HFOT)  / (1.  + ( 6 . ** . 5 /5 . ) * ALFA) 

ZU  = P(N2+1)  / AR( 1 ) 

ZD  = (HCRD  - (HRP(l)  - ZU) ) * 0.5 

ZX  = HR( 1)  - HFOT  - ZU 

C 

IF  (QPLM  .LT.  0.01)  THEN 
CSP  = CSU 
DTP  = TDF 
NSP  = NSU 
GO  TO  100 
END  IF 

CALL  PLUME 2 (CSU , QPLM , MSDT , NSDT , NSU , TDF , TL, ZI , ZJO , 

. BP, CSP, DTP, NSP, UP) 

100  CONTINUE 

C 

IF  (UP  .LT.  0.01)  THEN 
CSJ  = CSU 
DTJ  = TDF 
NSJ  = NSU 
GO  TO  200 
END  IF 

CALL  CLGJET  (BP , CSP , CSU , DTP , NSU , NSP , RDCT , TDF , TL, TW , UP , 

. CSJ  , DTJ  , NSJ  , HJ  , VJ ) 

200  CONTINUE 
C 

WRITE (7,700)  TIME , TU , DTP , DTJ , CSU , CSP , CSJ , NSU , NSP , NSJ , UP , VJ 
C 

RETURN 

700  FORMAT (1P12E10. 3) 

END 

C 
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SUBROUTINE  PLUME 2 (CSU , Q , QS ,QN,NSU, TDF , TL, ZI , ZJO , 

. BB , CSM , DTM , NSM , UM) 

C 

COMMON /SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 

. RDCT , RIO , SP , ST , TO 
C 

REAL  KK , LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL 
REAL  NSM, NSU 
C 

DIMENSION  WRK( 3 ) 

EXTERNAL  FPLM 
C 

DZ  =0.1 

TU  = TL  + TDF 

ZF  = 0.08  * Q**(2./5.) 

C 

FI  = Q*GRV*(LMD2+1. ) / (PI*GRO*GCP*TO) 

FS1  = QS*(LMD2+1 . ) / PI 
FN1  = QN*(LMD2+1.)  / PI 
C 

IF  (ZI  .LT.  ZF)  GO  TO  200 
C 

Z = ZI 

IF  (ZI  . GT . ZJO)  Z = ZJO 
C 

V = (9./5.*Fl*ALFA)**(l./3.)  * Z**(2./3.) 

W = 6 . /5 . * ALFA  * (9./5.*Fl*ALFA)**(l./3.)  * z**(5./3.) 

BB  = W / V 

UM  = (V*V)  / W 

DTM  = F1/(LMD2*GRV*W)  * TL 

CSM  = FS1  / (LMD2*W) 

NSM  = FN1  / (LMD2*W) 

IF  (Z  .GT.  ZJO-DZ/2 . ) RETURN 
C 

F2  = LMD2*GRV*W  * (DTM  - (LMD2+1 . )/LMD2  * TDF)  / TL 
FS2  = LMD2*W  * (CSM  - (LMD2+1 . )/LMD2  * CSU) 

FN2  = LMD2*W  * (NSM  - (LMD2+1 . )/LMD2  * NSU) 

C 

IF  (F2  .LT.  1.0E-04)  THEN 
BB  =0.0 
CSM  = CSU 
DTM  = TU 
NSM  = NSU 
UM  =0.0 
RETURN 
END  IF 
C 

WRK(l)  = W 
WRK(2)  = V 
WRK(3)  = F2 
C 

100  CALL  RKCAL(FPLM , Z , DZ , 3 , WRK) 

C 
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Z = Z + DZ 

IF  (Z  .LE.  ZJO-DZ/2 . ) GO  TO  100 
C 

W = WRK(l) 

V = WRK(2) 

GO  TO  300 

C 

200  F2  = FI 
FS2  = FS1 
FN2  = FN1 

DLTZ  = ((Fl/F2)**.5  - 1.)  * ZI 
Z2  = ZJO  + DLTZ 

V = (9./5.*F2*ALFA)**(l./3.)  * Z2**(2./3.) 

W = 6./5.  * ALFA  * (9./5.*F2*ALFA)**(l./3.)  * Z2**(5./3.) 
300  BB  = W / V 

UM  = (V*V)  / W 

DTM  = MAX(F2/(LMD2*GRV*W)*TU+TDF,  F1/(LMD2*GRV*W)*TL) 

GSM  = MAX(FS2/(LMD2*W)+CSU,  FS1/(LMD2*W) ) 

NSM  = MAX(FN2/(LMD2*W)+NSU,  FN1/(LMD2*W) ) 

C 

RETURN 

END 

C 
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SUBROUTINE  FPLM  (F,  X,  Y) 

C 

COMMON/ SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 
. RDCT , RIO , SP , ST , TO 
C 

DIMENSION  F(3) ,Y(3) 

C 

F(l)  = 2.  * ALFA  * Y(2) 

F ( 2 ) = Y( 3 ) * Y(l)  / Y(2)**3 
F(3)  = 0.0 
C 

RETURN 

END 

C 
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SUBROUTINE  CLGJET  (BB , CSP , CSU , DTP , NSU , NSP , RMX, TDF , TL, TW , UM, 
. CSJ  , DTJ  , NSJ  , HH , VM) 

C 

COMMON/SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 

. RDCT , RIO , SP , ST , TO 
C 

REAL  KK , LMD2 , MGAS , MGOUT , MSDT , MSTTL , NSDT , NSTTL 
REAL  NS A , NSAO , NSJ , NSP , NSU 
C 

DIMENSION  WRK( 3) 

EXTERNAL  FJET 


C 


C 


C 

C 


C 


C 


C 


C 


DR  =0.1 
TA  = TL  + TDF 


CSAO  = LMD2 

/ 

(LMD2+1 . ) 

* 

(CSP  - 

CSU) 

DTAO  = LMD2 

/ 

(LMD2+1 . ) 

* 

(DTP  - 

TDF) 

HAO  = BB  / 

6. 

** . 5 

NSAO  = LMD2 

/ 

(LMD2+1. ) 

* 

(NSP  - 

NSU) 

VAO  = UM  / 2 . ** . 5 

RIO  = GRV*(DTAO/TA)*HAO  / (VAO*VAO) 
RO  = 3 . ** . 5 * BB 

CSA  = CSAO 
DTA  = DTAO 
HA  = HAQ 
NSA  = NSAO 
R = RO 
RI  = RIO 
SGM  =0.0 
VA  = VAO 


IF  (RO  .GT.  RMX-DR/2 . ) GO  TO  200 

WRK(l)  = HAO 
WRK(2)  = RIO 
WRK( 3 ) = 0.0 

100  CALL  RKCAL  ( FJET, R, DR, 3 ,WRK) 

R = R + DR 

IF  (R  .LE.  RMX-DR/2.)  GO  TO  100 


HA  = WRK(l) 
RI  = WRK( 2 ) 
SGM  = WRK( 3 ) 


200  DMY  = EXP(-KK*ST/3.*SGM) 

VA  = VAO  * (RI/RIO  * R/R0)**(-l./3.)  * DMY 
CSA  = CSAO  * RO*HAO*VAO  / (R*HA*VA) 

DTA  = DTAO  * HAO/HA  * (RO/R*RO/R*RI/RIO)**( 1 . /3 . ) * DMY* DMY 
NSA  = NSAO  * RO*HAO*VAO  / (R*HA*VA) 


CSJ  = CSA  * ( (LMD2+1 . ) /LMD2 ) ** . 5 + CSU 
DTJ  = DTA  * ( (LMD2+1 . ) /LMD2)** . 5 + TDF 
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NSJ  = NSA  * ((LMD2+1.)/LMD2)**.5  + NSU 
HH  = HA  / (PI/2. )**.5 
VM  = VA  * 2 . ** . 5 
C 

RETURN 

END 

C 
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SUBROUTINE  FJET  (F,  X,  Y) 

COMMON/SMK/ALFA , BETA , DLTT , EJO , FW , GCP , GRV , GRO , HCRD , KK , 

. LMD2 , MG AS , MGOUT , MSDT , MSTTL , NSDT , NSTTL , PI , PMAS , PR , RCOAG , 

. RDCT , RIO , SP , ST , TO 

REAL  KK , LMD2 , MG AS , MGOUT , MSDT , MSTTL , NSDT , NSTTL 
DIMENSION  F( 3 ) ,Y(3) 

EJ  = EJO  * EXP (BETA* (RIO -Y( 2) ) ) 

F(l)  = (FW/2.  + (2.-  SP*Y(2)/2. )*EJ  - Y(l)/X  - SP*Y(2)*KK*ST/2 . ) 
/ (1.  - SP*Y(2)) 

F(2)  = (FW/2.  + (1.+  SP*Y(2)/2.)  * (EJ  - KK*ST/3 . ) 

- (1.+  2 . *SP*Y (2 ) ) * Y(1)/(3.*X)) 

/ ((1.-  SP*Y(2) ) * Y(1)/(3.*Y(2))) 

F(3)  = 1./  Y( 1) 

RETURN 

END 
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